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Abstract 

A deterministic mathematical model for the polymerization of hyperbranched 
molecules accounting for substitution, cyclization, and shielding effect has been 
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Introduction 

In this paper the mathematical modelling of hyperbranched polymers based 
on monomers of the AB 2 type is taken at hand, allowing the calculation of var- 
ious microstructural properties, using a numerical-mathematical scheme based 
on a graph-theoretical description of branching topology. . In previous work, 
scalar quantities like average molecular weight, polydispersity, degree of branch- 
ing, or gel point have been successfully predicted. These results were mainly 
achieved either by statistical theories of Stockmayer [TJ[2] and Flory [2j or, later, 
by stochastic simulations [4] and a generating functions approach El [7] . Most 
theoretical models that describe the polymerization of AB 2 monomer are based 
on the assumption that all functional groups of the same type are equally re- 
active, and react independently of one another. In other words, all B groups 
in the AB 2 monomer units in the polymer have the same reactivity. Yet, the 
equal reactivity is questioned if a substitution effect is taken into account [BJ. 

As regards the AB 2 problem with substitution, significant progress has been 
achieved in studies by Galina[5j [8], ChengJB] and Zhou[?l [9]. Here, the be- 
haviour in time of the degree of branching^, 6 , the dendritic and linear units[6] 
and the fraction of cyclized molecules [5] has been obtained. In the work by 
Galina[5 and Cheng [BJ the average microstructural properties have been com- 
puted using the generating functions approach. More recently, the 2-dimensional 
chain- lcngth/degree-of-branching distribution was obtained in the form of an 
analytical expression by Zhou[7] under the assumption that no cyclization re- 
action takes place. The intramolecular reaction of cyclization has indeed been 
considered in [5]. In view of its interesting gelation behaviour, the AB$ sys- 
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tern has been investigated using Monte Carlo simulation by Somvarsky et al. 
|10j Rather than substitution they were interested in the effect of shielding of 
monomer groups, leading to a shift of the gel point to higher conversion. They 
deal with the steric excluded volume effect by introducing so-called exponent 
kernels, thereby mimicking non equal accessibility of reactive functional groups. 
In the present paper we show that our newly developed numerical method can 
be easily adapted to capture the shielding effect in AB 2 system employing an 
exponent kernel. 

The new computational techniques that we will present in this paper are an- 
other contribution to the area of growing interest that is formed by the search 
after multidimensional distributive properties, where we previously treated poly- 
mer modification [TT] [TJ] and crosslinking polymerization [T31 [TTJ [T2] . In the 
realm of this trend here we will show how multidimensional distributive prop- 
erties can be recovered for the AB 2 polymerization system with substitution, 
cyclization, and shielding. The following reaction mechanism is considered: 



where a; is a number of terminal units, and y is a number of liner units in 
the hyperbranched polymer. The first two lines correspond to the polymeriza- 
tion reaction for substituted and non substituted groups. The third and forth 
lines denote the intramolecular reaction of cyclization. The fifth and sixth lines 
denote the reaction between cyclized and non-cyclized molecules. The rate co- 
efficients are not constant, and depend on x, y, chemical substitution ratio p, 
cyclisation factor A and scaling constant K t . Note, that we do not yet explicitly 
refer to shielding in this formulation. 
The paper is organised as follows. 
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Firstly, we discuss the main concepts of graph theory that has been used to 
describe molecular topologies of hyperbranched AB2 polymer. Graph theory 
is not only an important tool to build the mathematical model that is going 
to be solved. It also helps us understanding the connections between different 
distributive properties and, in addition, it allows expressing one property in 
terms of the other at the post-processing stage. Secondly, we show how the 
rules of topology evolution as induced by the chemical reactions, may be trans- 
formed into mathematical balance equations. The non-linear intcgro-differential 
equations will be discretized and solved numerically employing linearisation and 
projection methods. 

Thirdly, we demonstrate how the 'raw' numerical result, a three-dimensional 
distribution, may be post-processed in order to retrieve different distributive 
properties. The average properties and the chain length distribution obtained 
from the model, are shown to perfectly agree with findings of earlier investiga- 
tions [51 [7] for the case without cyclization. Furthermore, we show the first 
ever full multidimensional distributive properties for systems that experience 
both substitution and cyclization simultaneously. 

Fourthly, we show how the designed numerical method can be adapted to cap- 
ture the steric excluded volume effect by using the exponent kernels - the first 
time this concept is realized in a deterministic manner. 

1. Graph representation for hyperbranched polymers 

The topology of a hyperbranched polymer is now described in terms of graph 
theory. According to this approach, each monomer in the polymer of the AB2 
type has one node associated with it, while the chemical bonds between the 
monomers are represented by edges. A single polymer molecule is represented 
by a connected graph G with nodes of maximum degree three and at most one 
cycle. Moreover, the structure of a non-cyclized polymer molecule does not 
contain cycles at all, and is therefore called a rooted binary tree. In the present 
case of AB2 a cyclized polymer graph has exactly one cycle, see Figure [T] Now, 
it is useful to distinguish monomers according to their position in the topology 
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graph. Table [T] shows the correspondence between monomer position, graph 
theory terms, and commonly used short names 0[6]. 



Monomer 


Graph theory term 


Short name 


two reacted B groups, and unre- 


root node 


root, R 


acted A 






two unreactcd B groups, A re- 


node of degree one 


terminal unit, T 


acted 






one reacted and one unreacted 


node of degree two 


linear unit, L 


group B, A reacted 






all A and B groups reacted 


node of degree three 


dendritic unit, D 


polymer with no unreacted A 


graph with a cycle 


cycle, C 


groups 







Table 1: A graph theoretical frame of reference of AB2 polymer molecules. 




Figure 1: The graph on the left has a tree topology with root R, while the graph on the right 
has a cycle. The length of the cycle is equal to 4. 

The graph representation brings a lot of advantages along, for instance, using 
basic equalities from graph theory |14| one immediately infers the connection 
between the parameters given in Table [l] take place. In fact, it suffices knowing 
only two of the parameters in order to reconstruct the others, see Figure [2] 

Indeed, the dependence of the total number of units on the number of ter- 
minal and linear units is given by the rather trivial expression, 

N = 2T+L + C-l. (1.1) 

Now, the AB 2 polymerization process leading to branched topologies, can 
be viewed as a specific random graph process [15] that starts with n nodes and 
no edges, and subsequently at each step connects certain nodes according to the 
reactions taking place. 
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(Terminal, Total) > (Total, Linear) 



(Dendritic, Terminal) < (Linear, Dendritic) 

Figure 2: By knowing two of four parameters it is possible to recover the other two. 

2. Reaction mechanisms and mathematical model 

Let G = {G} be a collection of graphs (sometimes referred to as & forest |14| ) 
with, possibly repeating, components of the form G. Hence, G is a multiset. 
Employing this notation, we refer to G as to a single polymer molecule, while 
s € G represents all molecules instantaneously present in the system, and G all 
possible molecules that may be generated with the AE>2 polymerization process. 

Following the conventional approach for AB2 systems with substitution, two 
polymerization reactionsjU [6] and a cyclization[5] reaction are considered. Re- 
calling x, y denotes amount of terminal and linear units, the interpretation of 
the reactions in terms of the random graph process emerges as: 

• two arbitrary tree-components of G join together by connecting the root 
R from one component and a terminal unit T from another at rate K t , 

Px,y ~t~ T* x t yi y Px+x' — 1 .y+y' + l ; 

• two arbitrary tree-components of G (at least one of two is not strictly 
binary) join together by connecting the root R from one tree and a linear 
unit L from another at rate pK t , 

^x.y ~t~ T* x ' y' y Px+x' .y+y' — 1 ; 

• an arbitrary tree component of G receives a cycle by connecting the root 
R with a terminal unit T at rate XK t or by connecting root R and linear 
unit L at rate pK t , 
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• two arbitrary tree-components of G, of which one possesses a cycle, join 
together by connecting a root R from one molecule and a terminal unit T 
from another at rate K t ; or join by connecting root the R from one tree 
and a linear unit L from another at rate pK t , 

Px,y ~t" C x '.y' ^ ^x+x' — l.y+y' + l > 
"i" G x ' ,y' ^ Px+x' ,y+y' — 1 • 

Suppose we have a system of polymers s. Let each graph g E G has a number 
< F s (g) < 1 associated with it that corresponds to the relative frequency of 
occurrence of graph g within the system s £ G. In other words, for each selected 
topology g, F s (g) tells how probable it is finding g within a particular system s. 
As g contains full information on molecular topology, it is important to note that 
the whole system s, consisting of a large number of molecules, would have to be 
represented by an ensemble of the above-described graphs, which would require 
a huge amount of data to store. Thus, we will rather distinguish topologies only 
by three parameters (T,L,C), mapping G — > fl x {0, 1}, O = [0,oo) 2 . Here Q 
denotes the domain of all possible values for parameters (T,L). On the other 
hand, we expect s to be exclusively dependent on time by considering a function 
s(t). This reduces the original measure F s ( t ) to 

f c (x,y,t)eW, ce{0,l}, (2.1) 

which denotes the relative frequency of occurrence of a connected component 
with a number of terminal units between x and x + dx, a number of linear units 
between y and y + dy, and c cycles in the polymerization system at time t. 
We refer to W = L 2 (il) x C[l,oo), as a short way to describe mathematical 
properties of the distribution f(x,y,t). Namely, L 2 (Q) suggests f c (x,y 7 t) is 
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integrable over the first two variables x,y; C[l,oo] mark it is continuous in 
time, t. We will deliberately use this shortening in the paper. For instance, 
when referring to a mapping W — > W, a, law that sets correspondence between 
two distributions from W; a mapping W 2 — > W - correspondence between a pair 
of distributions and another distribution from W ; and W — > R - correspondence 
between a distribution for W and a real number. 

Passing to the limit dx, dy — > the evolution of f c (x, y, t) with time may be 
described with a system of differential equations. Two equations are derived, 
for cyclized polymers c = 1, and for non-cyclized ones, c = 0, 

df c (x,y,t) 



at 

subject to the initial conditions 



= (C c f)(x,y,t), c = 0,l (2.2) 



f (x,y,0) = 5(x-l)6(y), 
fi(x,y,0)=0, 



(2.3) 



where S is a delta function and operator C :W — > W implements the effect 
of the reactions listed at the beginning of this paragraph on the distribution 
f c . The operator fully determines the dynamics of the system and, as it will 
be shown further on, has a non-linear structure. Now, we start with a set of 
basic operators, typical for polymerization, that also act on f c and eventually 
will serve as building blocks to define a more complicated structure of C 

Definition 2.1. The shift operator T Xg yg g W — > W takes a function f e W 
to its translation T Xs ^ Vs f , 

T Xs , ys f{x,y,t) = f(x-x s ,y-y s ,t), (x s , y s ) E O. (2.4) 



We normally associate the operator (2.4) with propagation-like reaction mech 



anisms, where at each firing of the reaction one node is added to the connected 
component G. 

Example 2.1. Let us consider a simple propagation reaction Pi iV + M 
Pi, y+ i, meaning that the lengths of the linear chains Pi iy are growing at constant 
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rate K t . The corresponding differential equation has the form 
§j.f( x > y> *) = K t (To,if{x, V, t) - /(at, y, t) ^ 



The differential equation from Example 2.1 may be viewed as a partial case of a 
chemical master equation |16j . Note besides that any chemical master equation 



may be represented as a linear combination of shift operators (2.4). 

Definition 2.2. Integral operators fj,, fj, xi [i y : W — > K, take a function f to its 
partial moments /z/, /z^/, fi y f: 

(m/)(*) = / f(x,y,t)dxdy, 
n 

iVxf)(t) = J xf(x,y,t)dxdy, (2.5) 
n 

= J y.f( x ,y,t)dxdy, 

a 

Example 2.2. Let f c {x,y,t) € W be a known evolution of the topology for 
an AB 2 polymerization system as described in Section [7} Then, as a direct 
consequence of the equality the total mass time profile for the non-cyclized 

polymers (c = 0) is given by 

m(t) = (2 fa + pL y - (i)fo(x, y, t). 

Definition 2.3. Convolution is an operation that takes two arguments f,g^zW 
to a product f * g defined as 

x V 

f(x, y, t) * g(x, y,t) = J J /(£, r), t)g(x - y - r], t) d^dn (2.6) 
o o 

It may easily be checked that the convolution is a symmetrical f * g = g * f , 
and a bilinear a ■ f * g = a ■ (f * g). atl operation. Note, that the operation 
exhibits special behaviour regarding the 0-moment: 

H(J * g) = Hf > M (2.7) 

The convolution is an important tool when modelling step-growth polymeriza- 
tion [11] . Firing of such a reaction in terms of a random graph process means 
that any two arbitrary disconnected components from G become connected by 
a new edge. This is illustrated by the following example. 
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Example 2.3. Let us consider a polymerization reaction described by the reac- 
tion equation P x .y + Px'.y' Px+x' .y+y' , where the reaction rate is not depen- 



dent on x,y,x',y'. In contrast to to Example (2.1) there are multiple choices 
to select the reactants, so P XjV is produced for a fixed combination x, y. The 
corresponding differential equation has is formulated as, 

J^/fo y> l ) = Kt (/fo V' f ) * /fo y> f ) " /fo V' *) ■ ^/fo y> *)) • 



The differential equation presented in Example (2.3) can be also viewed as a 
Smoluchowski coagulation equation [17] . 



Now, we are ready to define C(f c (x,y,t)), c = 0, 1, as appearing in (2.2). 
In analogy to the ideas expressed in Examples 1 2 . 1 1| 2 . 3] we assemble the operator 

A)(/c) = ^* ^ T o,i /o * a;/o - (Mx/o + A*/o x) ■ f - Xxfo^j + 

pK t (t-l-x f * yf - (nyf + pf y) ■ f - Xyf ^j ; (2.8) 

£i(/ c ) = KtiTo^xfx *f - xfxfifo) + pKtiTi-tyft * f - yfiftfo) + XK t (xf + 

(2.9) 

Here as before, the subscript denotes the balances for cyclized 1 and non-cyclized 



molecules. The equation (2.8) consists of two similar lines that correspond to 
R + T and R + L reactions. Here, p > is the chemical substitution factor. As 



similar to Example 2.3 To,i/o *xfo denotes the production term. Here, however, 
the weight x has been used as the reaction rate is proportional to the amount of 
terminal units. The shift operator T^i was employed, since with every reaction 
firing, one extra linear unit appears. The consumption term — (Mx/o+m/o x ) ' fo 
shows that each polymer might participate in the reaction by contributing either 
an R unit at rate pfo ■ xfo or a T unit at rate p x fo ' fo- Finally, the cyclization 
R + T occurs at rate XK t xfo, where A > is a constant factor. The second line 



of (2.8) may be analyzed in an analogous manner with as a major difference 
that it denotes the substituted reaction R + L this time. 



It may be seen that the balance equation for components without cycles (2.8 1 
is presented in a closed form, so the differential equation for fo may be solved 
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independently of f\. In the case /o is known, /i may be resolved by integration 

9fi(x,y,t) . 

— Q~ t — = £i{fc{x,y,t))- (2.10) 

For this reason we subsequently consider two differential equations instead of 



the full system (2.2). The numerical approach is similar for the cases and will 
first be introduced for Cq. 

2.1. Numerical treatment 

We start with the time discretization. In view of the non-linear nature of 



(2.2), it is important to employ an implicit time integration scheme, also known 



as Rothe's method. Let f(x^y,t k ) be an approximation to /q, the solution of 



(2.2). Then, an implicit first order time integration scheme is given by 

f(x, y, ife+i) = f(x, y, t k ) + T k £ f(x, y, t k+1 ), k = 0, 1, 2, . . . 
or in a more compact way, 

{T k £o-I)-f(x,y,t k+l ) + f(x,y,t k ) = 0, k = 0,1,2,... (2.11) 
where / is an identity operator and r k = t k +i — t k is a time step. Now, there are 



two main issues in solving the equation (2.11) for f(x,y,t k+ i). On each time 
step we deal with: 



• Linearization: find a linear equation that mimics the behaviour of (2.11) 
locally with respect to (a;, y). 

• Approximation: find an approximate solution that satisfies the linear equa- 
tion up to certain error tolerance. 

Let £ / : W 2 — > W be the Frechet derivative [18] of Cq at point f{x, y, t k ), 

C Q fh = K t I T ,i (/ *xh + h* xf ) - fi x f ■ h - {p,f + X)x ■ h I + 

V ' x (2-12) 

P K t \ Tx-i (f*yh + h* yf) ~n y f-h- (fif + \)y ■ h\ . 

Here, the direction h G W of Frechet derivative should not be confused with 



the time step r k . Now, a functional Newton's method may be applied to (2.11 ) 
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that defines a sequence {/ s }o° converging to the root f(x,y,tk+i) of equation 



(2.111 



f 



s+1 



r 



(2.13) 



Here, we assume that the time step index k is fixed and the initial estimation 
f°{x, y, tk+i) is given. Now, with each iteration s = 0, 1, 2, . . . , the new estimate 
f s+1 is obtained as an upgrade of f s with a certain correction term h s that solves 
the equation 

(T k C' Qfs - I)h s = (r k C - I)f s +f(x,y,t k ). (2.14) 

Although the equations ( 2.13| 2.14) contain Co, they are linear with respect to 
the correction term h s . However, the estimated solution f s+1 is expressed only 



implicitly as shown in (2.14). Thus, we proceed by approximately solving (2.14) 



employing a collocation projection into a Gaussian basis|llj. 

Let a(h, v) : W 2 — > K, b(v) : W — > M be a bilinear and a linear form, 



respectively, being associated with the right and left hand sides of (2.14) 



a(h, v) :=< {T k C f , - I)h, v >, 
b(v) :=< ( Tk C f s -I) + f tk) ,v>, 



(2.15) 



where < ■, ■ > is an inner product 



< f(x, y),g(x, y) >= J f(x, y)g{x, y) dx dy. 
u 



(2.16) 



The equations (2.13 2.14) may be reformulated in a weak sense 



(2.17) 



f s+1 = f s - h s , h s € W; 
a(h,v) = b(v), Vw € W. 

The idea behind a 'weak' formulation is the following: instead of requiring the 



residual in the original equation (2.14) to be zero, we rather demand a zero 



inner product of the residual and all test functions v £ W . This is formulated 



in the second line of (2.17) 



Finally, in order to perform the transition from the dimcnsionally infinite W 
to the finite degrees of freedom (DoF), we construct a system of n basis function 
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centres (xi,yi) € f2 and connectivity parameters a.^. Two n-dimensional bases 
in W are defined as follows 



(f>i(x,y) = e 



-cr Xli (x-Xi) -iy y ,i(y-yi) 



a. a > 0, i = 1,2,. 



ipi(x,y) = 5(x - Xi)S(y - y*), i = 1,2, . . . ,n. 
W n = span{(f>i,<p2, • ■ • , <M c W, 
W tes i = spani^x, i/j 2) ipn} C W, 



(2.18) 

(2.19) 
(2.20) 
(2.21) 



The basis functions (2.18) are used for expansion of the approximation h, f s € 



W n , where h approximates h and, to the weak solution h G W, 

n 
n 



(2.22) 



(2.23) 



while (2.191 is a source of test functions v used in (2.171. Thus, we arrive at a 



discrete collocation scheme for equation (2.14) 



a(h n ,ip j )=b(Tp j ), j = l,2,...,n (2.24) 
For the sake of brevity we write a = (ctj), (3 = (pi) referring to the column 



vectors of coefficients as defined in (2.22 2.22) and employ this matrix notation 



henceforth. Thus, for instance, relation (2.24) represents a system of n linear 



equations and by substituting (2.22) and (2.15) into (2.24) one may bring it to 
the following matrix form 

M°a = b (2.25) 

where a is a column of unknown coefficients defining the approximation to the 
correction term h s , whileM is a square n x n matrix, 



T k K t i f 0) i (Cpf x + C fx(3 ) - (t x (3I n - (A/3 + X)f x 



TkpKAT^lCpTy + C. 



fi y (3I n - (A/3 + \)T y - L, 



(2.26) 
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Here, /„ is an identity matrix of size n, 



T — 4 _1 A 



(2.27) 



is a discrete approximation to the operator (2.4) 



<f>j(xi - x s ,yi - y s ) 



(2.28) 



(2.29) 



Matrices 

T x = A _1 diag{xi,a;2, . . -,x n }A, 

T y = -A -1 diag{j/i,2/2, ■ ■ -,y n }A 
represent the multiplication with weights x or y, respectively. Matrices Cp, Cf B 
represent the convolution with the known approximation f s or its weighted form 



wf s associated to the s-th iteration of the Newton process (2.13) on the fc-th 
time step. 



A~ l C, 



( c )i,j = J2 Mj(%i,yi)*(t>k(xi,yi). 



(2.30) 



Analogously to (2.26), by expanding (2.15) the column vector at the right 



hand-side of (2.25) may be obtained as a matrix expression 



b = ( 2>n,Kt{T 0tl CpT x - fi x f3I n - (pp - X)T X 



Finally, the partial moments fi^, (t) , fly (t) of the approximated distribution may 
be obtained by the relation 



<t>i(x,y)dxdy, 

'n 

p*(t) = q T T%(3. 



(2.32) 
(2.33) 



The considerations expressed so far show how the approximation to the 
distribution fo{x,y,t) of cyclized molecules may be retrieved by solving the 
differential equation involving Cq. The situation regarding f\{x,y,t) is even 
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simpler as the differential equation (2.10) has a linear form, provided Jo is 
known. Thus, a linear matrix transform of the coefficient column vector a tk on 



time step is sufficient to obtain data on time step t^+i ■ 

a t k+1 = OL tk + T k Ml k+i a tk+l 



(2.34) 



Ml = 2K t (C p f x - q T (3f x ) + pK t {C p f y ~ q T (3T y ) + \K t (T x (3 + P T y f3) (2.35) 



Until now, only a first order time integration has been considered. The 
numerical tools derived in the previous section may provide sufficient data for 
higher order time integration techniques. For instance, the computational codes 
belonging to multistep methods for stiff systems require subroutines computing 



the coefficients vector for the discretized right hand side of (2.2), a c = C c f3, C c : 



R n , and its Jacobian matrix J a (3. Here, we strongly benefit from the 



specific construction of the test function space (2.19 2.21 ) that satisfies 



< f{x,y),ipj{x,y) >= f(xj,yj), ipj € W test - 
Furthermore, the explicit formulas for the coefficients are obtained analogously 



to the expressions (2.26 2.31) 



£o/3 = ( 2Kt { T 0A CpT x - (j, x I n - fiT x j + 

+ K' t (fl-lCpTy - fJ,yI n - pTy^j(3\ 

j£ o (3 = 2K t ^T 0i i (Cpf x + Cf^^j - [i x I n - fif s 



(2.36) 



In the case of C\ the corresponding differential equation (2.10) is linear, hence 

U/3 = J to P = Ml, (2.37) 

In conclusion, at this stage we have developed a fully discrete set of formulas that 
may be implemented into computational code in a straightforward manner or 
treated employing existing stiff numerical integrators (e.g. MATLAB odel3s). 
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2.2. The numerical algorithm 

In order to provide a summary of the instructions listed in the previous 
section, the numerical algorithm is depicted in the diagram of Appendix 1. The 
discretization techniques presented fully describe the approximated systems by a 
collection of column vectors f3 tk , k = 0, 1, 2, . . . that correspond to each time step 
tk- Regarding space approximations, the basis function centres are enumerated 
with a single integer number, so notwithstanding the dimensionality of two, the 
problem could be parametrized using a vector of n coefficients at a single time 
point. Furthermore, Newton's iteration process is established at every time step 
tk, so we introduce an upper index s to refer to the sequence of estimations in 
the Newton's process: j3 8 t . 

Before the algorithm starts, one must set the basis parameters Xj, yi, o~ x y . 
This could be realized either on the basis of the previous simulations in or- 
der to refine the results in certain parts of the domain, or as a logarithmi- 
cally distributed system that covers the whole domain of possible values (x,y), 
[0,x max ] x [1, y max \. There is no simple way to compute optimal values for 
o~x,y, although the parameters should be dependent on the distance between 
two adjustment basis functions. Since the collocation approach tends to give 
the smallest absolute error at the basis function centres, the approximation er- 
ror could be a posteriori evaluated locally, by repeating the simulations using 
additional checkpoint basis functions in between the original ones. This princi- 
ple is illustrated in Figure [3] The the approach of residual subsamphngpjj] is 
applied to refine the mesh. 

• If the local error is greater than a certain tolerance Tol ac id, the interme- 
diate checkpoint node A is accepted in the new system. 

• By default all original basis centres are accepted in the new system. 

• If the local error in four intermediate nodes is smaller than a predefined 
tolerance Tol remove , the original node R situated in between is not ac- 
cepted into the new system. 
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Figure 3: The principle of adaptive mesh refining is based on an error evaluation at interme- 
diate nodes x. The original nodes system is denoted by o 



The matrix expression derived in the paper could become computationally 



intensive. Nevertheless, as it is shown in (2.26) and (2.31), the matrices M and 



b, dependent on time and a certain intermediate result, are to be recomputed on 



every step and consist of constant numbers (2.27 2.28 2.29). that only depend 



on the fixed basis functions parameters, and hence should only be computed 
once at the initialization of the algorithm. 



3. Results and post-analysis 

The numerical scheme described in the previous paragraph shows many ad- 
vantages when comparing with existing simulation techniques. It provides a 
representation of the full distribution, that is accurate enough to perform data 
mining even in the regions with very small values. Despite high precision, the 
method remains computationally inexpensive. As a result we are able to extract 
very detailed morphology related properties by post-processing. For instance 
information on cycle length, that was typically associated with Monte Carlo 
simulations before, is obtained in a fully deterministic manner. 

By following the numerical scheme we generate full time profiles for non- 
cyclized and cyclized molecules fo(x, y, t), fi(x, y, t) that correspond to the model 
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012345678 
Terminal units, log 10 



Figure 4: The level-lines of weighted frequency distribution xyfo(x,y,t finc i) as obtained from 
simulations. Different values for the substitution factor are considered p S {0.1,1,10}. For 
each case 8 level-lines are plotted io-( 4 - 5 +°- 5fc ) , k = 1, . . . , 8. 




Figure 5: The level-lines of weighted frequency distribution xyfo(x,y,t flrll i) as obtained 
from simulations. Different values for the cyclization parameter are considered A £ 
{10" 3 ,10^ 4 , 10~ 6 }. In each case 10 level-lines are plotted 10 _(4 - 2+ 5 fe) , k = 1, ... ,8. 
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Terminal units, log 10 



Figure 6: The level-lines of cyclized molecules frequency weighted distribution xyfi (x, y, t enl i) 
as retrieved from fo(x,y,t). Different values for the substitution factor are considered p £ 



{0.1, 1, 10}. In each case 7 level-lines are plotted 10 



-(5+ife) 



,10. 



9 
A 

Conversion of A groups, c 

Degree of freedom for approximation 

%i i Hi 

Tolerance for Newton iterations e 
Initial time step 

computational time, single core CPU: 
for quadratures (to be done once) 

for time integration (to be done for each parameter setup) 

Table 2: Parameters used in simulations 



0.1, 1, 10 

0, IO" 6 , IO" 5 , 10" 4 , 10~ 3 

0.9, 0.99, 0.997, 0.998, 0.999, 0.9999 

1200 

{1, 2, 3, 2 fc , fc = 2 . . . 8Zo.g 2 10} 
{2.7(xi - Xi-iy 2 ,i = 1, . . . ,n} 

io- 14 
io- 10 

35 min. 
5 min. 



(2.2) up to a conversion of A groups c = 99.95%. The full set of parameters is 



given in Table [2j The effect of two input parameters is studied in detail: 

• p the ratio between the reactivity of substituted and non-substituted 
groups, see Figure |4j 

• A the ratio of the cyclization rate to the polymerization rate, see Figure [5j 

The distribution for cyclized structures fi(x,y,t) is retrieved by integrating 
( 2~10| , Figure [6] 
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The results shown in Figures 4[6 already provide an idea how the substi- 
tution factor p influences the topologies. Indeed, varying the values for p is 
observed shifting the two dimensional distribution diagonally along the diag- 
onal line logs + logy = const. This equally holds for non-cyclized molecules 
Figure |4j and cyclized molecules Figure [6] The cyclization factor A shifts the 
distribution along log x — logy = const as shown in Figure [5] This - at first sight 
uncomplicated - behaviour of the two dimensional distribution nevertheless has 
dramatic implications for the scalar and distributive properties that are to be 
inferred by post processing. 

In view of high amount of statistical properties that are extractable from 
the simulation results f c , we group them into four categories: time dependent 
scalars and 1,2,3-dimensional distributions. 

3.1. Time dependent scalars 



The time dependent scalars are obtained from the moments (2.5) of the 



density f c (x,y,t) as a convenient consequence of basic graph properties (1.1). 
All time dependent values will be presented with respect to conversion c G [0,1) 
of the free A groups, so K t coefficient may be chosen arbitrary. Indeed, the 
time dynamics for the conversion of A groups c(t) ; the number of free A groups 
n A (t); the number of linear rii(t), terminal or dendritic njy(t) units; the 

fraction of cycles nc(t) or the degree of branching may be expressed as follows: 

C W = v(fo(x,y,t) + fi(x,y,t)); 

n A (t) =l-c(t); 

n L (t) = p y fo(x,y,t); 

n T (t) = p x f (x,y,t); (3.1) 

n D (t) = n T (t) - n A (t); 

n c(t) = pfi(x,y,t)/c(t); 

db(t) = — . 

\ / n D +0.5n L 

Here we illustrate the most important of the listed properties. The Frey's de- 
gree of branching[20] db(t) turns out being strongly affected by the substitution 
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factor p, while the cyclization parameter A does not influence the results signif- 
icantly, Figure [7] On the other hand, the fraction of cyclized molecules nc(t) 
is dependant on A but not on p, Figure [8] This results are in good agreement 
with previous findings [SJ |5] . 




0.2 0.4 0.6 0.8 1 

Conversion of A groups 



Figure 7: Degree of branching db(t) as a function of conversion of A groups c(t) in polymer- 
ization of AB2 monomers reacting with a substitution effect for B-groupsp. 
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Figure 8: Fractions of cycle-containing molecules nc(t) versus conversion c(t) as calculated 
from the simulation. The parameter A controls the extent of cyclization. The substitution 
effect p does not affect this property. 
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3.2. 1- dimensional distributions 

One-dimensional distributions are obtained by evaluation of line integrals 
of fc{x,y,t en d) for a fixed time point t en d- Distributions of molecular weight, 
degree of branching, or cycle lengths may be obtained in this fashion by defining 
a collection of lines that keeps a certain property constant, see Figure [9j 
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Terminal units, log^ 

b ) 




012345678 
Terminal units, log 10 

c) 

Figure 9: Three collections of lines that keep constant values of a) chain length , b) de- 
gree of branching, c) length of a cycle. The level lines of weighted frequency distribution 
x yfo{ x ^y^end) with A = 0, p = {0.1, 1, 10 arc plotted in the background for reference. 

For instance, in order to recover the amount of molecules with a certain 
length n from fo(x, y, t en d), one has to consider all combinations of points (x, y) 
that satisfy 2x + y — 1 = n (see Figure [9^,), and sum the corresponding values 
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of fo- Mathematically, this is equivalent to evaluation of a line integral along 
the collection of lines ld n (y) = 0.5(n — y + 1), y € [0, n — 1], 

n-l 



ld(n) = 



1 



fo{ld n (y),y,t end )dy 



(3.2) 



The molecular weight distribution n 2 ld(n) for various substitution factors and 
no cyclization is shown in Figure |10| The effect of different conversion values as 



the distribution is approaching its asymptotic values is shown in Figure f 1 The 
effect of approaching a straight line in a double logarithmic plot as a system 
approach the gelation point was predicted before for a simpler svstems|21j. It is 
also important to consider the cases without cyclization, as an analytical solu- 
tion is known [7], which may be employed to validate the numerically obtained 
results. As shown in Figures \lO\l 1| perfect agreement is observed over a span of 
16 orders of magnitude. 

The effect of the cyclization factor A on the chain length distribution for a 



fixed conversion value and substitution factor is shown in Figure 15 




10 

Chain length 



Figure 10: Normalised double weighted chain length distribution n 2 ld(n). The predicted 
values are indicated by the dashed lines and compared to the analytical solution, the solid 
lines. The effect of three levels of chemical substitution p S {0.1,1,10} is illustrated under 
the assumption that no cyclization takes place. 
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Chain length 



Figure 11: Chain length distribution ld(n) as predicted by the numerical method, indicated by 
the solid lines, as compared to the analytical solution, in dashed lines. The effect of different 
values for conversion c S {0.9, 0.99, 0.998, 0.999} is illustrated under the assumption that no 
cyclization takes place A = 0. The asymptotic value for conversion 1 is illustrated by the 
dash-dot line. 



The length distributions of polymers with different degree of branching b € 
(0,1] is found by integrating over the collection of lines defined by lb(x) — 
2^^-x, x £ [l,oo) (see Figure ^p) s Analogously to (3.2) we obtain 

oo 

bd(b) = J(2x + y)f c (x, l b (x),t end ) dx. (3.3) 
o 

Note, we intentionally use the chain length 2x + y as a weight while integrating 
in order to magnify the contribution of longer molecules, that have much more 
monodisperse branching distribution. Although, the weight was used we still 



observe prolonged tails in the branching distribution bd(b) Figure 13 This fact 
has to be accounted for when one uses average an value instead of a distribution 
to describe the branched topology of the AB2 system. We also observe that the 
cyclized molecules contribute less to the tails than non-cyclized molecules. 

In order to obtain information concerning cycle lengths, we employ an idea 
formulated by to Panholzer et al.[22]. It was shown that the expectation of 
a distance depthi^ n from the root to an i-labelled terminal unit is statistically 
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Figure 12: Chain length distribution ld(n) as predicted by the numerical method. The solid 
lines denote cases with a cyclization factor A = 10 -4 , for different values of conversion. The 
dashed line depicts the non-cyclization case for reference. The effect of different values of 
conversion c e {0.99, 0.997, 0.999, 0.9999} is illustrated. 



connected to the number of terminal units a; in a strictly binary tree: 

2(2» + l)(2a;-2 t + l) Q ( 2x x ~_f ) 
depth z . x = — j^r l,i = 0,...x. (3.4) 

x + z UJ 

Now, by proportionally extending the length of each path by a number of linear 
units and averaging over the label index i, an expected cycle length cl(x, y, Ct ype ) 
is obtained, 



c/(x, 7/, C-type) 



x + 0.5y - 1 



^2depth ijX . 



(3.5) 



The level lines of (3.5 1 are illustrated in Figure [9J:. Note that it is important to 
know whether the reaction that caused cyclization added an edge between the 
root and a terminal unit Ct ype = or between the root and a linear unit Ct ype = 



1. As follows from the model (2.8), at each instant of time X(xfo + pyfo) new 
cyclized molecules appear. The instantaneous cycle length distribution of these 



molecules may be extracted according to principle (3.5). In order to retrieve the 



total accumulated distribution of cycle lengths we have to additionally integrate 
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0.2 0.4 0.6 0.8 1 

Branching degree 

Figure 13: Normalised Frey's degree of branching distribution as obtained from the simulations 
for two levels of chemical substitution p £ {0.1, 10}. The overall distribution (the dashed lines) 
is plotted together with its two components, the distributions for cyclized and non-cyclized 
polymers (solid lines). The rate of cycle formation is A = 10 — 4 . 



over the whole time span [0,t e nd]. This requires the evaluation of the double 
integral, 



tend 



cd(n) = J J Xxf (cl(x,y),T)d(x,y)dT+ 

cl{x,y,0)— n 



^pyfQ(cl{x,y),T)d{x 7 y)dT. (3.6) 



cl(x ,y,l)—n 



As in previous cases we evaluate the results (3.6) in a numerical manner. The 



result is depicted in Figure [l4| That it turns out to be possible to retrieve the 
cycle length distribution proves the great potential of this deterministic method, 
especially when taking the few dimensions considered into account. It provides 
sufficient information to obtain morphology related properties that are usually 
only attainable with Monte Carlo simulations. 
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Figure 14: Normalised cycle-length distribution as obtained from the integration of 
fc(%, y, tend)- The effect of different values for cyclization factor A £ {10 — 8 , 10 — 5 , 10 -4 } 
is illustrated. 




Figure 15: The effect of the different rates of cycle formation A 6 {10 6 , 10 5 , 10 4 } on the 
chain length distribution. The dashed line depicts the case without cycle formation. 

3.3. 2- dimensional distributions 

The distribution f c ix^y,t en( i) denotes the relative frequency of structures 



with x terminal units, y linear units, and c £ {0, 1}, as depicted in Figures 4[6 



It is interesting to change the coordinate system of terminal units versus linear 
units (x, y) to another system. For instance, the coordinate system chain length 
versus degree of branching (d, n) may provide information that would remain 
unnoticed otherwise. In this case the following transformation of coordinates is 
required, 

x — y — - — 

x+Q.5y ^ 7) 

y — » 2.x + y — 1 



According to (3.1) the first line of the coordinates transform replaces x with 
Frey's degree of branching, while the second line replaces y with chain length. 
The resulting distribution is depicted in Figure |16| One may observe that 
the distributions are centred around the average property: Frey's degree of 
branching. However, they remain fairly broad for chain lengths smaller than 
10 3 , while accompanied by extremely narrow peaks for longer chain lengths. 
This fact remains hidden when only the average degree branching is considered. 
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Figure 16: Two-dimensional distributions of chain length and degree of branching. Three 
values for the substitution factor are considered, p £ {0.1, 1, 10}. In each case 8 level lines are 

plotted, l0~( 3+ 2 k \ k = 1, . . . , 8. The dotted line depicts a scalar degree of branching for the 
whole system. 

3.4- Dynamic evolution of 2D distributions 

As an example of 3-dimensional distribution obtained from the simulations, 



we present a time trajectory for distributions shown in Figure 16 Here f c {x, y, t) 



is not only considered at the end time point t en d but similar to ( 3.1 ) on the whole 
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time interval t £ [0, t en( j\ as depicted in Figure 17 Note, as in the previous 
subsection we use double-weighted chain length in order to highlight the long 
molecules present at extremely low concentrations. 



p = 0.1 




Figure 17: Molecular structures with different degree of branches and molecular weight distri- 
bution at 10 time-slices. Two values for the substitution factor are considered, p £ {0.1, 10}. 
In each case 10 level lines are plotted, lQ-( 2 +°- 5k ) , fc = 1, . . . , 10. 
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4. Non-linear substitution as a consequence of shielding 

As has been pointed out by Somvarsky et. al. [TO] H3] the reactivity of 
a functional group may be dependant on the position of the group within the 
molecule's topology due to a steric excluded volume effect. The amount of non- 
shielded terminal x' and linear y' units may be defined by employing a power 
law 

x' oc x u \ 

(4.1) 

¥ oc y u , 

where u) is the constant that describes the shape of the molecules, as a measure 
for the reactive surface having a fractal character [TU]. The qualitative investi- 
gation of the kernels with fractional and integer values for w has been performed 
by Wattis [21] ■ The exactly solvable cases correspond only to integer values for 
oj [2~T]. In case of a fractional to, the method of generating functions [2"4] would be 
particularly difficult to apply as the transformation of the kernel into a domain 
of generating functions would imply a fractional differentiation, whose definition 
usually requires an integral transformation, which in this case would return the 
original (unsolved) problem. In the previous sections, we have already consid- 
ered the case of possibly substituted, but equally accessible groups oj = 1 and 
as we shall see further on that there are no particular difficulties to expand the 
approach to cases of fractional lo. 

In the current paragraph the opposite marginal case of spheric particles 
oj = | is studied. The numerical scheme allows accounting for this shielding 
effect by applying only minor alterations to the original equations. It suffices 



reformulating the weight operators (2.291 in the following manner 



f£ =A- 1 di ag {yf,y^,...,y^}A 

By taking analogous steps as before we arrive at an approximate solution to 
the polymerization problem with shielded groups. The chain length distribution 
at different levels of conversion is presented in Figure [18] As a result, assuming 
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non equal accessibility for reactive functional groups leads to a slowing down of 
reaction rates, and consequently shorter molecules are formed. 




Chain length 

Figure 18: The chain length distributions at different levels of conversion, c S 0.9, 0.99, 0.999 
as obtained by numerical simulation accounting for the shielding effect, w = 2/3 (depicted with 
dashed line). The chain length distributions for identical parameters, but without shielding 
u = 1 are depicted with solid lines for reference. 

Conclusions 

The novel numerical method capable of capturing two-dimensional popula- 
tion balances rigorously, was developed as a result of clear mathematical rea- 
soning. In contrary to prior results, the method reveals a full two-dimensional 
distribution. The method represents distributions as a linear combination of 
Gaussian basis functions, and finds their expansion coefficients, finally resulting 
in an approximate solution of the population balance equation. The accuracy 
of more than 15 orders-of-magnitude was observed when comparing with the 
analytical results. 

The distributive properties of hyperbranched ABi system with a substitution 
has been obtained by applying the method to a population balance model. 
In contrast to prior modelling studies of AB 2 polymerization, the distributive 
properties are obtained in multiple dimensions and maintain high resolution, 
giving access to information of large molecules that are only present in extremely 
low concentrations. In addition, an intramolecular cyclization reaction has been 
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accounted for in the model. The influence of the chemical substitution and 
shielding as well as various cyclization rates on the molecular topologies has 
been studied. 

The numerical results obtained for the case of no cyclization turned out to be 
in excellent agreement with the analytical solution present in the literature [7]. 
For the cyclization case we observe agreement with averaged properties obtained 
by GalinaJS]. However, the newly developed method is able to reveal much 
more detailed information. Interestingly, the simulations show that the degree 
of branching is strongly dependent on the substitution factor, while the chain 
length is predominately determined by cyclization rate. Low cyclization rates 
cause small amounts of molecules with long cycle length and vice versa, high 
cyclization causes more molecules to have cycles, but shorter in length. It was 
also found that the branching distribution has a broad non-trivial shape that 
is non-symmetrical for cases with substitution effect. The distribution can be 
relatively well described by Frey's degree of branching index (FDB) for long 
molecules where it deviates less from the expected value. In contrast, for short 
molecules the deviation increases and FDB only poorly describes topologies. It 
has also been observed that the degree of branching for molecules with cycles is 
generally more narrowly distributed than for non-cyclized molecules. 

Cyclization turns out to have two different effects on the chain length dis- 
tribution. Higher rates of cyclization tend to shorten the tail at very long chain 
length, but at the same time the chain length distribution decays more slowly 
in the mid range. The effect of the cyclization on chain length and branching 
distributions may be explained by the fact that cyclized molecules react at a 
lower rate, which attributes a kind of 'memory' to the polymerization system 
and blurs the distributions towards the earlier time stages. 

It was also shown how the chemical substitution implemented as a linear 
factor, can be adapted to model the shielding effect. The shielding effect is 
implemented by a non-linear kernel with a fractional power, and assumes only 
those functional groups are accessible for reaction that are situated on the sur- 
face with dimension of fractal character. 
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The developed algorithm provides a fast and efficient way for predicting 
multidimensional distributive properties for hyperbranched polymers. This was 
successfully demonstrated on the example of AB 2 type polymerization account- 
ing for the substitution effect and cyclization. Although the paper is more than 
just a study of a particular molecular system, we have presented a methodol- 
ogy that is capable of handling a wide spectrum of polymer related problems 
and is especially helpful in cases of more than one dimension and convolution 
equations. 
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Appendix 1 



k = k + 1 







Initialization 


• 


Fix the degree of 




freedom, n 


• 


Fix Xi,yi,<T Xt i,a y ,i 




i = l,2, ...,n. 


• 


Compute constant 




matrices, equations 




1.28-1.30, 


• 


set P to to corre- 




spond initial condi- 




tions (1.4) 





Resolving 

• Compute matrix 1.27 

. s = 0, fi" = fi tk 

• Compute column vec- 
tor 1.32 and solve linear 
system 1.26 

. =/3> -a 

• if T x T y \p°+ l - P»\ >€, 
set s = s + 1 and repeat 
last two steps; 
otherwise, /3 tk+1 = f) s+1 . 



refine b, 



Finalization 

• Error estimation. 

• If necessary increase 
DoF / adapt basis 
centers and restart the 
process. 

• Postprocessing and 
output. 
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